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James R. DeBonis 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

A computational fluid dynamics code that solves the compressible Navier-Stokes equa- 
tions was applied to the Taylor-Green vortex problem to examine the code’s ability to ac- 
curately simulate the vortex decay and subsequent turbulence. The code, WRLES (Wave 
Resolving Large-Eddy Simulation), uses explicit central-differencing to compute the spatial 
derivatives and explicit Low Dispersion Runge-Kutta methods for the temporal discretiza- 
tion. The flow was first studied and characterized using Bogey & Bailley’s 13-point dis- 
persion relation preserving (DRP) scheme. The kinetic energy dissipation rate, computed 
both directly and from the enstrophy field, vorticity contours, and the energy spectra are 
examined. Results are in excellent agreement with a reference solution obtained using a 
spectral method and provide insight into computations of turbulent flows. In addition the 
following studies were performed: a comparison of 4th-, 8th-, 12th- and DRP spatial dif- 
ferencing schemes, the effect of the solution filtering on the results, the effect of large-eddy 
simulation sub-grid scale models, and the effect of high-order discretization of the viscous 
terms. 


Nomenclature 


Q"n 

spatial finite difference coefficients 

bn 

spatial filter coefficients 

C, C* 

actual and numerical wave speed 

e t 

total energy 

p 

pressure 

Qi 

heat flux vector 

t 

time 

t* 

nondimensional time, t* = t/(L/Vo) 

Ui 

velocity vector 

U, V, w 

streamwise, transverse and spanwise velocity components 

X, y , 2 

cartesian coordinates 

D 

spatial finite difference operator 

Ek 

kinetic energy 

E (E k ) 

spectral density of kinetic energy 

L 

reference length 

M 

Mach number 

Re 

Reynolds number 

Sij 

strain rate tensor 

cd 

deviatoric strain rate tensor, Sf,- = Sij — SkkSij 

T 

temperature 

V 0 

initial velocity 

&rm Pm 

low dispersion Runge-Kutta coefficients 

Sij 

Kronecker delta 

At 

time step 
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At* 

e 

ei 

£3 

K 

P 

V 

LVi 

n 

p 

a 

T ij 

c 


nondimensional time step, At/(L/Vo) 
kinetic energy dissipation rate 

kinetic energy dissipation rate computed from deviatoric strain rate 

kinetic energy dissipation rate computed from pressure dilatation 

angular wave number 

kinematic viscosity 

spectroscopic wave number 

vorticity vector 

volume 

density 

filter damping coefficient 
stress tensor 
enstrophy 


subscripts 

0 initial value (t* = 0) 


superscripts 

spatial filtered 
Favre filtered 
sgs sub grid scale 


I. Introduction 

The understanding and prediction of turbulent fluid flow continues to be one of the pacing items in 
the physical sciences. The bulk of turbulent flow predictions are done using Reynolds-Averaged Navier- 
Stokes (RANS) methods. RANS solves a time-averaged form of the Navier-Stokes equations and the effect 
of the unsteady turbulent motion on the time-mean ffowfield is approximated using a turbulence model. 
RANS methods perform well for “well-behaved” flows where the models have been calibrated, e.g.. attached 
boundary layers and mildly separated flows. They fail in more complex flows beyond the bounds of the 
modeling assumptions, e.g.. large separations, highly anisotropic flows, shock boundary layer interactions, 
etc. RANS methods rely completely on the model for the prediction of the turbulent ffowfield. To reduce or 
remove the reliance on modeling, methods that resolve some or all of the turbulent motion are gaining usage. 
Direct numerical simulation (DNS) solves the unsteady Navier-Stokes equations and computes all scales of 
the turbulent motion and eliminates modeling entirely. Large-eddy simulation (LES) directly computes the 
large-scale turbulent structures, those that contain most of the turbulent kinetic energy, and uses a model to 
capture the effects of the smaller unresolved scales, minimizing the model’s influence. The success of both 
DNS and LES rely on the ability of the numerical scheme to accurately compute the unsteady turbulent 
motion. 

The Taylor-Green vortex (TGV) is a canonical problem in fluid dynamics developed to study vortex 
dynamics, turbulent transition, turbulent decay and the energy dissipation process. 1 The TGV problem 
contains several key physical processes in turbulence in a simple construct and therefore is an excellent case 
for the evaluation of turbulent flow simulation methodologies. The problem consists of a cubic volume of 
fluid that contains a smooth initial distribution of vorticity. Periodic boundary conditions are applied to all 
boundary surfaces. As time advances the vortices roll-up, stretch and interact, eventually breaking down into 
turbulence. Because there is no external forcing the small-scale turbulent motion will eventually dissipate 
all the energy in the fluid and it will come to rest. 

In this work we explore the TGV problem using the Wave Resolving Large-Eddy Simulation (WRLES) 
code 2,3 to gain a better understanding of the interplay between the numerical methods, grid resolution and 
turbulence, and to validate the code. An initial set of solutions was generated for the First International 
Workshop on High-Order Methods in Computational Fluid Dynamics, which was held at the American 
Institute of Aeronautics and Astronautics’s Aerospace Sciences Meeting in Nashville Tennessee on January 
7-8, 2012. The TGV problem was one of several cases that participants could solve and submit their results 
for comparison. Additional solutions have since been generated for this paper. 
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II. Code Description 


The code used in this study, WRLES (Wave Resolving Large-Eddy Simulation), is a special purpose large- 
eddy simulation code that uses high-resolution temporal and spatial discretization schemes to accurately 
simulate the convection of turbulent structures. 2,3 The code solves the compressible Favre-filtered Navier- 
Stokes equations on structured meshes using generalized curvilinear coordinates. It is written entirely in 
Fortran 90 and utilizes both Message Passing Interface (MPI) libraries 4 and OpenMP compiler directives 5 
for parallelization. 

A. Governing Equations 

The basis of large-eddy simulation is the separation of the large and small scale turbulent fluctuations. The 
large-scale turbulence is computed directly using the Navier-Stokes equations. The small scale turbulence is 
modeled. Since the large-scales carry the vast majority of the turbulent kinetic energy the technique offers 
promise for accurate turbulent simulations. Since the small scale turbulence serves to dissipate the large 
scales and is isotropic, a simple model can be constructed to impose the effect of the small scales on the rest 
of the flow. 

A filtering process is applied to the continuity, momentum, and energy equations. The resulting equa- 
tions are comprised of resolved and unresolved terms. The resolved terms in the filtered equations directly 
correspond, in form, to the unfiltered equations. The resulting LES expressions for conservation of mass, 
momentum, and energy are: 


dj_ 

dt 


d. {ni) = 0 


(i) 


8 ' 

dt 


(pe t ) + 


+ w, + p) = ay ( T « + T « ' ') 

d d 

— ( pUilt + Uip) = — [(ujTij + h ( r,/' ) - (q t 


\ 


(2) 

(3) 


The overbar, ( ) represents a spatially filtered quantity and the tilde, ( ) represents a Favre-averaged 
quantity; f = ( pf)/p 

The unclosed terms from the LES momentum and energy equations are the sub-grid scale stress, t‘? s , 
and the sub-grid scale heat flux, q sgs . In a traditional LES approach a sub-grid model would be used to 
compute and q S9S . Both the standard Smagorinsky model 6 and Lilly’s dynamic Smagorinsky model' are 
included in the code. Many practitioners forego sub-grid modeling, relying on the dissipation implicit in the 
numerical scheme to dissipate the energy at the small-scales. This approach, sometimes called implicit LES 
(ILES), is used extensively here. 


B. Numerical Method 

The Favre-filtered Navier-Stokes equations are discretized and solved in strong conservation form using the 
explicit numerical methods described below. 


1. Time Discretization 

Time discretization is performed using a family of explicit Runge-Kutta schemes written in a general M - stage 
2 -N storage formulation. 8 

dfm = eXmdfm— 1 T AtD (/m— 1) ( 4 a) 

fm = fm—1 + Pmdfm (4b) 

for m = 1 . . . M, and where the initial value at the start of the iteration is /o = /” and the final value is 
/m = f n+1 - The coefficient oq is typically set to zero for the algorithm to be self-starting. The operator D 
is the spatial finite difference operator. 

The order of accuracy and dispersion characteristics of the temporal discretization can be varied by 
changing the number of stages and the coefficients, a m and (3 m ■ The 4-stage, 3rd-order algorithm of Carpenter 
and Kennedy 9 was used for all computations in this study. 
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2. Spatial Discretization 

Spatial discretization is performed using central difference methods. The central difference stencil can be 
written for an arbitrary stencil size 


d_l 

dx 


N 1 

^ ' ~^^[ a n(fi+n ~ fi—n )] 


( 5 ) 


where the half width of the stencil, the number of points on either side of the central point, is N. The total 
number of points in the stencil is 2N + 1. The spatial discretization scheme can be varied by simply changing 
the width of the stencil and the coefficients, a n . Standard 4th-, 8th- and 12th-order central differences (St04, 
St08 & Stl2) and a 13-point dispersion relation preserving (DRP) scheme from Bogey & Bailly 10 (BB13) 
were used in this study. 

Central differencing is used for the spatial discretization because of its non-dissipative properties. Tur- 
bulent structures that are well resolved will be accurately convected without being damped. The standard 
schemes are designed to minimize the truncation error for a given stencil size. The DRP scheme is designed 
to minimize the dispersion error (phase error) and is formally 4th-order accurate. 


3. Filtering 

While the lack of dissipation in central differencing is beneficial for the accuracy of turbulent simulations, 
these schemes are not stable without some form of damping. To ensure a stable solution without adversely 
affecting the resolving properties of the scheme, solution filtering is used. Solution filtering is a low-pass 
filter that selectively adds dissipation to damp the high-wavenumber/unresolved structures that lead to 
instabilities but leaves the low- wavenumber/ well-resolved structures untouched. The filter must be properly 
matched to the differencing scheme, so that the filter removes only those waves that are not properly resolved. 

N 

fi = fi + V X! b ™fi+n (6) 

n——N 

For the standard central difference schemes, Kennedy and Carpenter’s 11 4th-, 8th- and 12th-order filters 
(KC04, KC08 & KC12) are used. For the 13-point DRP scheme the matching filter developed by Bogey and 
Bailly is used (BB13). 10 


4- Properties of the Numerical Scheme 

Fourier analysis can be used to evaluate the ability of numerical schemes to resolve wave motion. 12 The 
technique analyzes the numerical scheme on the one-dimensional convection equation and can provide the 
scheme’s behavior in terms of dissipative and dispersive errors. The effect of the temporal discretization is 
not considered here. 

The analysis shows that central difference schemes have no inherent dissipation, and thus are ideally 
suited for LES calculations. They do however produce dispersive errors. Figure la shows the phase error, 
written as the ratio of numerical phase speed to actual phase speed, c*/c, versus normalized angular wave 
number, (kAx), for the schemes used in this study. For poorly resolved waves, high values of k Ax, the 
scheme is unable to accurately resolve the wave. 

The damping function of the filter can also be expressed in terms of k Ax. Figure lb shows the damping 
behavior of the KC04, KC08, KC12 and BB13 filters. The damping is minimal until it reaches a “cutoff” 
wave number, where it rapidly increases, effectively removing all the structures of higher waves numbers. 
The strength of the filter can be reduced from the default behavior by adjusting a coefficient (ranging from 
0 to 1) multiplying the dissipative term. 


III. Problem Setup 

The TGV problem can be run using a variety of flow conditions and initial conditions. The conditions 
and post processing used here were specified by the organizers of the AIAA First International Workshop 
on High-Order Methods in Computational Fluid Dynamics. 
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(a) Error in phase speed for the finite difference stencils 



(b) Damping functions of the explicit filters 


Figure 1: Properties of the numerical scheme 


A. Flow Conditions 

The domain consists of a cube with sides of length 2nL. Within the domain an initial distribution of velocity 
and its corresponding vorticity is specified by the following relations. 


u 

v 

w 

V 


Vi) sin (y) COS fy) cos fy) 

— *«■(!) 

0 


Po + 


PqVq 

16 




( 7 ) 

( 8 ) 

(9) 

( 10 ) 


The Reynolds number selected for the workshop was 1600. Because WRLES solves the equations in 
dimensional form additional flow conditions, given in table 1, were selected to provide an incompressible flow 
at the given Reynolds number. 


Table 1: Flow conditions 

Quantity Value 

Reynolds number, Re 1600 

Mach Number, M 0.1 

Reference length, L 0.005 ft. 

Temperature, T 530 R 


B. Data Processing 

The primary method for evaluating the TGV solutions is examining the rate at which the fluid dissipates 
the kinetic energy in the domain. The integrated kinetic energy within the domain is computed by 
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The kinetic energy dissipation rate (KEDR) can then be computed by differencing Ej- in time. 


e(E k ) 


dEk 

dt 


(12) 


This KEDR, e, is written 
energy in the domain. 

A second “theoretical” 
£, within the domain 


as a function of E k to emphasize that it is computed directly from the kinetic 
kinetic energy dissipation rate can be computed from the integrated enstrophy, 


c 


1 

PqVL 



(13) 


It can be shown that for incompressible flow, the enstrophy is directly related to the kinetic energy dissipation 
rate through a constant. 


e (C) = 2— c 

Po 


(14) 


The difference between the directly computed KEDR, e(E k ), and the enstrophy based KEDR, e(£), will be 
examined. 

Two additional quantities were computed to test the incompressible assumption. For a compressible 
flow where the bulk viscosity can be taken to be zero, the “theoretical” dissipation rate is the sum of two 
quantities; the contribution to the KEDR, obtained from the deviatoric strain rate tensor 




and the contribution to the KEDR, obtained from the pressure dilatation term. 

£3 = 


1 

PoO 


L 


diii 

p- — ail 


dxj 


(15) 


(16) 


For low Mach number flows £3 should be negligible and the “theoretical” dissipation rate can be approximated 
using the enstrophy based KEDR (equations 13 & 14). 


C. Grids 

The grids were generated using a Fortran code that also computed the initial solution at t* = 0. The grids 
are equally spaced cartesian grids of 64 3 , 128 3 , 256 3 , and 512 3 points. In order to maintain the high-order 
of accuracy at the boundaries of the domain, 11 additional planes are added beyond the x = +7 rL, y = +7 tL 
and z = +nL boundaries and periodic conditions are enforced over a range of 6 points adjacent to the 
boundaries. This insures that all points within the domain are computed using the full finite difference 
stencil and filter stencil, and that the resolution of the scheme is maintained. 


IV. Results 

For the baseline and numerical scheme investigations, grid resolution studies were performed using grids 
of 64 3 , 128 3 , 256 3 , and 512 3 points. For each differencing scheme and grid, the coefficient that multiplies the 
effect of the filter was halved until a minimum value was found that provided a stable solution with minimal 
dissipation. It was found that this minimum filter coefficient was the same for all grid resolutions for a given 
scheme. The nondimensional time step used for the 64 3 case was At* = 3.385 • 1CP 3 . The time step was 
halved for each doubling of the grid dimensions. 

The 64 3 and 128 3 cases were run on a single processor desktop workstation with six cores. The larger 
cases were run on the NASA Pleaides high performance computing system using 8, 8-core processors for 
the 256 3 cases and 46, 8-core processors for the 512 3 cases. Approximate wall clock time for the runs is 
given in table 2. The WRLES code was originally written for the 13-point DRP scheme. The full 13-point 
stencil is always solved and there is no reduction in work for the low-order schemes. The other schemes are 
implemented by changing the differencing stencil coefficients and zeroes are used where necessary for the 
low-order schemes. For this reason there is no CPU time difference between schemes. 
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Table 2: Computing resources 


grid size 

time step (At*) 

machine 

cores 

wallclock time (hrs.) 

64 3 

3.385-10 -3 

desktop 

6 

0.5 

128 3 

1.693- 10 -3 

desktop 

6 

9 

256 3 

8.463-10 -4 

Pleiades 

64 

40 

512 3 

4.231T0 -4 

Pleaides 

368 

130 


Numerous solutions were run to support the following studies 

• Baseline study using the 13-point Bogey and Bailley DRP scheme (BB13) 

• Comparison of numerical schemes: standard 4th (ST04), 8th (ST08) and 12th (ST12) order central 
differencing and 13pt. DRP scheme (BB13) 

• Examination of the energy spectra 

• Effect of the filter 

• Effect of sub-grid stress model 

• Effect of the order of accuracy of the viscous term derivatives 

A. Baseline 

A baseline set of simulations was performed using the BB13 scheme. Four grid levels were run and a filter 
coefficient of 0.05 was used for all the cases. Iso-contours of the z-component of vorticity (figure 2) illustrate 
the evolution of the flow as described by Brachet in reference 13. At the earliest times the flow behaves 
inviscidly as the vortices begin to evolve and roll-up. Near t* = 7 the smooth vortical structures begin 
to undergo changes in their structure and around t* = 9 the coherent structures breakdown. Beyond this 
breakdown, the flow is fully turbulent and the structures slowly decay until the flow comes to rest. 

The change in the kinetic energy over time is shown for all four grid levels in figure 3. Little difference is 
seen between the grid levels in this plot; although a close-up of 10 < t* < 20 indicates that the coarser the 
grid, the less energy it contains as time evolves. 

Figure 4 shows the turbulent decay process. The directly computed KEDR, e(Ek), is given in figure 4a. 
Reasonable agreement with the reference solution 14 is show for all grid levels. The largest discrepancies are 
at the peak dissipation rates, 7 < t* < 12. At t* = 9, where the dissipation rate peaks, the 128 3 , 256 3 , and 
512 3 grids are all within two percent of the reference solution (the 64 3 grid has an error of 9.8 percent). 

The KEDR computed using the enstrophy integral, e(£), is given in figure 4b. In theory the enstrophy 
based KEDR should be equivalent to the directly computed KEDR, but it is clear that there is a large 
discrepancy in the peak dissipation rate for the lower grid levels (46 percent for the 64 3 case) that improves 
with grid resolution. The 256 3 (within 2.5 percent) and 512 3 (within one percent) grids are in excellent 
agreement with the reference solution. 

The contributions to the KEDR based on compressible flow assumptions are shown in figures 4c and 
4d. The primary contribution to the KEDR is through ei and the data shows that it is identical to the 
incompressible e(£) predictions. The pressure - dilatation contribution, € 3 , is two-orders of magnitude smaller 
than ei and can be neglected. The €3 curve for the 256 3 case differs from the other cases beyond t* > 8. 
This behavior is unexplained and does not affect the results. These results confirm the incompressible flow 
assumption, and it will be used for the remaining analysis and discussion. 

Contours of the vorticity norm on the constant x-plane, x = — 7 rL, at time t* = 8 are shown in figure 
5. On the 64 3 grid, the vorticity is smeared over a large area, the structure is jagged and not well defined, 
and the peak levels are low relative to the higher resolution grids. As resolution is increased, the structure 
becomes less smeared, smoother and more defined, and the peak vorticity level increases. There are minimal 
differences between the 256 3 and the 512 3 solution and those plots closely resemble the reference solution 
shown by van Rees. 14 The increased resolution of the vortical structures with grid resolution seen here, 
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Figure 2: Iso-surfaces of z-vorticity from the BB13 scheme on the 128 3 grid 




(a) complete simulation, 0 < t* < 20 


(b) close up of 10 < t* < 20 


Figure 3: Evolution of kinetic energy 
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(a) directly computed KEDR, e(-Efc) 



(c) deviatoric strain tensor based KEDR, ei 



(b) enstrophy based KEDR, e(£) 



(d) pressure dilatation based KEDR, £3 


Figure 4: Kinetic energy dissipation rate; comparison of grid resolutions using the BB13 scheme 
directly corresponds to the improved prediction of KEDR computed from the enstrophy (equation 13) shown 
in figure 4b; enstrophy is simply the square of the vorticity magnitude. 
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(a) 64 3 


(b) 128 3 


(c) 256 3 


(d) 512 3 


1 3 5 7 ^9 11 13 15 | 

Figure 5: Contours of the vorticity norm |cj| pr, at t* = 8 on the plane x = — nL, in the region y = 0 to ^ 
and z = ^ to nL; comparison of grid resolutions using the BB13 scheme 


The discrepancy between e(Ek) and e(£) provides valuable insight into turbulent simulations. Despite 
the huge variation in resolution (a factor of 512 in the number of grid points between the 64 3 and 512 3 
simulations) the kinetic energy in the domain (figure 3) and the directly computed KEDR (e(£ , j c )) (figure 
4a) are well predicted at all grid levels. However, at lower grid resolutions, the KEDR computed from the 
enstrophy (e(C)) is significantly lower than the directly computed KEDR. In other words the dissipation 
due to the vorticity (i.e. turbulence) does not account for all the dissipation present in the simulation. 
An additional mechanism is present and it plays a significant role at low resolutions. This mechanism is 
numerical dissipation. 

For the WRLES code the numerical dissipation is generated solely by the spatial filter. As the grid is 
refined the resolution of the vortical structures increases, increasing e(£) and the effect of the filter decreases, 
reducing numerical dissipation. The fact that e(Ek) is well predicted at all grid levels indicates: 1) the two 
sources of dissipation in the calculation are in proper proportion to each other and 2) the filter correctly 
mimics the physical dissipation process of the unresolved scales. This validates the use of these numerical 
methods for implicit large-eddy simulations. The addition of an explicit sub-grid model will be examined in 
a later section. 

B. Comparison of Schemes 

A comparison of the four spatial differencing schemes on all grid levels is shown in figure 6. The directly 
computed and theoretical KEDR are shown. The results are consistent with the baseline study (figure 4). 
The directly computed kinetic energy dissipation rate agrees reasonably well for all grid levels and schemes, 
with the exception of the 4th-order scheme on the 64 3 grid, where the combination of lower-order scheme and 
coarse grid can not resolve the flow adequately and the resulting KEDR curve exhibits a different character 
from the others. The enstrophy based KEDR is under-predicted for all schemes until the 256 3 grid is reached. 
The lower order schemes all approach the correct e(£) levels with grid refinement, but the 4th-order scheme 
requires additional refinement beyond 256 3 . The ST12 and BB13 schemes, both utilize 13 point stencils and 
their predictions are very similar with some differences occurring near the peak dissipation and beyond. 

C. Examination of Spectra 

Kinetic energy spectra were computed to show the evolution and decay of the turbulent structures. Spectra 
were computed along each grid line and averaged. A saw-tooth shaped odd-even grid point oscillation in 
the spectra was seen for all cases. Brachet 13 also reported this behavior and averaged the spectra over twice 
the band width ( 2nLAv = 2) to remove the oscillation. This averaging was used here as well. 
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(a) 64 3 grid 


(b) 64 3 grid 




(c) 128 3 grid 


(d) 128 3 grid 


Figure 6: Kinetic energy dissipation rate; comparison of finite differencing schemes. Left column, directly 
computed KEDR, e(Ek), right column, enstrophy based KEDR, e(£) 
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(e) 256 3 grid 


(f) 256 3 grid 




(g) 512 3 grid 


(h) 512 3 grid 


Figure 6: (continued) Kinetic energy dissipation rate; comparison of finite differencing schemes. Left column, 
directly computed KEDR, e(Ek), right column, enstroplry based KEDR, e(£) 
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The progression of the energy spectra with time on the 512 3 grid is shown in figure 7. At t* = 2 all 
the energy is confined to small wave numbers. As time progresses, the energy begins to cascade down to 
smaller and smaller scales. The energy in the smallest scales peaks around t* = 10. This corresponds to the 
maximum energy dissipation rate. Beyond this time, the energy decays for all wave numbers. 



Figure 7: Evolution of the energy spectra for the BB13 scheme on the 512 3 grid 

The spectra for t* = 12 for both the BB13 and ST12 schemes on all grid levels are shown in figure 8. For 
the lower grid resolutions one can see that the energy near the end of the spectra “tails off,” falling below 
the higher resolution curves. This reduction in energy is the result of the filter dissipating the energy at the 
highest wave numbers in the solution, where the differencing scheme is producing dispersion errors. Refining 
the grid moves this effect to higher and higher wave numbers. The difference between the ST12 and BB13 
schemes is also illustrated by the spectra. The DRP schemes trade formal order of accuracy for improved 
resolution of waves. Careful comparison of the close-up plot shows that for a given grid resolution the energy 
in the ST12 solution begins to fall away from the higher resolution curves sooner than the BB13 solution. 
However the difference is slight. 



(a) energy spectra at t* = 12 (b) close-up of energy spectra at t* = 12 


Figure 8: Energy spectra at t* = 12 for the BB13 and ST12 schemes 
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D. Effect of the Filter 


Artificial/numerical dissipation is present in some form in all CFD codes. This dissipation is necessary to 
obtain stable solutions, but too much dissipation can have an adverse effect on the results. This is especially 
true for DNS and LES analyses where the numerical dissipation can easily damp the turbulent structures 
that are critical for an accurate solution. For the WRLES code, artificial dissipation is incorporated through 
the filter. There are two key parameters that control the filtering process: 1) the shape of the filter response 
(unique to each filtering scheme) and 2) the damping coefficient, a. For this study the BB13 scheme was 
used to solve the TGV problem on a 128 3 grid. Three filtering schemes were examined: 1) the baseline case 
using the matching BB13 filter and the minimum stable damping coefficient of 0.05, 2) the BB13 filter and 
the maximum damping coefficient of 1.0 and 3) the 6th-order Kennedy & Carpenter filter, KC06, and a 
damping coefficient of 0.05. 

Figure 9 shows the evolution of the kinetic energy for the 3 filtering schemes. Increasing the damping 
coefficient maintains the shape of the filter response and simply increases the amount of damping across all 
wave numbers. This is reflected in the fact that the shape of the kinetic energy curve is maintained, but 
the amount of energy is reduced. This is also illustrated by examining the kinetic energy spectra at t* = 12 
(figure 10). The shape of the two BB13 filter curves is very similar, but the a = 1.0 curve contains less 
energy across the entire spectra. Changing the filter to KC06, reduces the cutoff wave number; this will 
remove larger turbulent structures. By removing larger structures the solution initially contains less energy 
(5 < t* < 10). But because the large scale turbulent structures are responsible for dissipating energy the 
turbulent dissipation is lessened at later times and the kinetic energy is greater than the baseline solution 
beyond t* = 15. The spectra for the KC06 curve (figure 10) shows similar energy content at the low 
wavenumbers, but the energy contained at the higher wave numbers is reduced. 



Figure 9: Evolution of kinetic energy, BB13 scheme on the 128 3 grid; effect of filter order 


The kinetic energy dissipation rates are shown in figure 11. These plots confirm the above observations. 
For both modified dissipation cases, the additional damping causes an increase in KEDR for the first half of 
the simulation. Near the peak dissipation, the KC06 case shows significantly less e(Ek) because there are less 
large scale structures. The added dissipation across all wavenumbers for the maximum coefficient case also 
shows a reduction in e(E k ), but not to the same extent. As expected the enstrophy based KEDR indicates 
that both cases significantly increase the amount of numerical dissipation in the simulation; the KC06 case 
being the most dissipative. 

In section C the pronounced reduction in the magnitude of the spectra at the highest resolved wave 
numbers was attributed to the filter. This effect is examined more closely here. The error in the spectra on 
a given grid was obtained by comparing that spectra to the spectra on the highest resolution grid (512 3 ), 
|E(/c)5i2 — E(fc) n |. This error, for the BB13 scheme on the 256 3 grid, is shown in figure 12. The error is plotted 
against the normalized angular wave number. The damping function of the filter is plotted for comparison. 
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Figure 10: Energy spectra, BB13, 128 3 grid; effect of filter parameters 




(a) directly computed KEDR, e(-Efc) (b) enstrophy based KEDR, c(C) 

Figure 11: Evolution of kinetic energy and kinetic energy dissipation rate; effect of filter parameters 
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For t* < 4 (figure 12a) there is no discernible trend; perhaps because the flow is still transitioning and the 
spectra are not fully developed. At later times (figure 12b) a clear trend is seen: 1) at low wave numbers 
the error is low, 2) at approximately kAx = 2 the error begins to rise steadily, and 3) the error mimics the 
shape of the filter’s damping function. Clearly, the filter is removing the energy beyond the nominal cut-off 
wave number of kAx « 2 while leaving the well resolved wave numbers untouched. The t* = 20 curve does 
not follow these trends and the reason is unknown. The observations made here are also valid for the other 
grid resolutions. 



(a) t* < 4 (b) t* >5 

Figure 12: Error in the energy spectra for the BB13 scheme on the 256 3 grid compared to the damping 
function of the filter 


E. Effect of sub-grid model 

All the solutions presented so far have been without any explicit modeling of the sub-grid scale turbulence. 
Large-eddy simulation sub-grid scale models were examined using the BB13 scheme and the 128 3 grid. Sub- 
grid models mimic the effect of the small scale turbulence by generating additional dissipation that is intended 
to affect the resolved turbulent structures. The standard Smagorinsky model and the dynamic Smagorinsky 
model are compared to the no model case in figure 13. Examining the KEDR curves(figures 13a & 13b), both 
sub-grid models perform similarly. Early in the simulation the flow contains large coherent structures and a 
fully turbulent flow has not yet developed. There is no small scale turbulence and the sub-grid model is not 
needed. However, the sub-grid models do not recognize this and produce dissipation that adds to the overall 
dissipation in the simulation. This is evidenced by the increased e(E k ) for t* < 8. Beyond this point, the 
flow is fully turbulent and the sub-grid model’s dissipation acts appropriately, damping the large turbulent 
structures that are the primary dissipation mechanism. This effect is seen in the reduction of e(E k ) from 
the no sub-grid model case for 8 < t* < 12. At later times the models perform similarly to the no model 
case. Because the models add dissipation to the simulation and damp the large scale turbulent structures, 
they significantly reduce the enstrophy based KEDR (figures 13c & 13d), and increase the disparity between 
e(E k ) and e(Q. 

The dynamic model produces a jagged e(E k ) curve at the peak dissipation rates and also produces slightly 
higher levels of e(£) at that location. The dynamic model locally varies the Smagorinsky coefficient based 
on the resolved flowfield. In general, this feature results in a less dissipative sub-grid model and also allows 
for the backscatter of energy; energy flowing from the sub-grid scales, back to the resolved scales. The 
jagged appearance of the curve is due to this backscatter effect. This was confirmed with a simulation where 
backscatter was prevented and the curve became smooth. 

The effect of the sub-grid model on the resolved structures can be visualized using vorticity contours 
(figure 14). The no model case has the highest level of vorticity and a slightly jagged appearance to the 
contours. The sub-grid model solutions maintain the general characteristics, but the peak vorticity levels 
are reduced and the contours have been smoothed. 
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(a) directly computed KEDR, e( Ef.) (b) close-up of directly computed KEDR, e(Ej.) 




(c) enstrophy based KEDR, e(£) (d) close-up of enstrophy based KEDR, e(£) 


Figure 13: Kinetic energy dissipation rate; comparison of sub-grid scale models. 
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1 3 5 7 9 11 1315 

Figure 14: Contours of the vorticity norm M^, at t* = 8 on the plane x = —i tL, in the region y = 0 to 
and z = ^ to nL; comparison of sub-grid models using the BB13 scheme on the 128 3 grid 


F. Effect of high-order viscous terms 

It is standard practice in high-order Navier-Stokes codes to limit the high-order schemes to the convective 
terms and compute the viscous terms with lower order methods. In the case of WRLES and many other 
codes, fourth-order differencing is used for the viscous derivatives. It has been long assumed that the accuracy 
of the convective terms is critical for resolving the turbulence and the accuracy of viscous terms is much less 
important. This treatment reduces computational cost and complexity. 

It was suggested at the High-Order CFD workshop that the author investigate the effect of using the 
same order of accuracy for both convective and viscous derivatives. In this investigation, the standard 
12th-order scheme was used for the convective terms and solutions were obtained with both 4th-order and 
12th-order accurate viscous terms. Figure 15 shows the results. e(E^) is only slightly effected by the high- 
order viscous terms and no clear trend is visible. However, figure 15b clearly shows that the higher order 



(a) directly computed KEDR, e(-Efc) (b) enstrophy based KEDR, e(£) 

Figure 15: Kinetic energy dissipation rate; effect of high-order viscous terms 
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viscous terms reduce the e(£). The lower-order representation under-predicts the magnitude of the derivatives 
used in computing the stress tensor. Using the more accurate higher-order representation results in larger 
magnitudes for and hence a more dissipative term in the equations. 

V. Conclusions 

A numerical study of the Taylor-Green vortex problem was undertaken using central-differencing based 
finite difference methods. The compressible Navier-Stokes equations were solved using both standard central 
differencing and Dispersion Relation Preserving differencing schemes. Grids ranging from 64 3 to 512 3 were 
used. 

For the baseline method, Bogey & Bailley’s 13-point Dispersion Relation Preserving (DRP) scheme, the 
kinetic energy dissipation rate (KEDR), computed directly from the kinetic energy in the domain, e{E^) was 
predicted reasonably well at all grid levels, with the largest discrepancies occurring near the peak dissipation 
rates. A second “theoretical” KEDR based on the enstrophy in the domain, e(() was significantly under- 
predicted on the lower resolution grids. Increasing grid resolution significantly improved the results. Results 
at the finer grid levels (256 3 and 512 3 ) were in very good agreement with a reference solution generated with 
a spectral method. 14 

For many solutions the directly computed KEDR is well predicted and the enstrophy based KEDR 
is under-predicted. This indicates that: 1) the turbulent structures are under-resolved, 2) the numerical 
dissipation from the filter plays a significant role in energy dissipation, and 3) the current numerical approach 
is well suited for implicit large-eddy simulation. 

Several spatial differencing schemes were compared: the 13-point DRP scheme of Bogey & Bailley and 
the 4th, 8th and 12th order central differencing schemes. Results were as expected with the higher-order 
schemes achieving grid convergence faster than lower-order schemes. To accurately resolve the details of 
the turbulent structures, high-order schemes are recommended. However, lower-order schemes were shown 
to adequately resolve the directly computed KEDR. The DRP scheme and 12th order scheme had very 
comparable solutions; this is expected as they have the same size stencils. The DRP scheme slightly extends 
the resolved range of the turbulent energy spectra, but this difference is small. 

Solution filtering was used to maintain numerical stability and the minimum stable value of the damping 
coefficient was used for each solution. The filter was shown to behave as designed, damping the poorly 
resolved high wavenumber structures and leaving low wavenumbers untouched. A pronounced “tailing off” 
of the energy spectra at the highest wave numbers was shown to be due to the filter. Increasing the filter’s 
damping coefficient increased the KEDR. Increasing damping by using a lower-order filter with a smaller 
cut-off wavenumber was shown to actually reduce the KEDR because some of the large-scale turbulent 
structures responsible for the dissipation have been removed. 

The effect of large-eddy simulation based sub-grid models were examined at the under-resolved 128 3 grid 
level. Both the Smagorinsky and dynamic Smagorinsky models were compared. The sub-grid models were 
found to: 1) add too much dissipation where the flow is not fully turbulent, 2) improve the solution near the 
peak dissipation rates, 3) dissipate the resolved turbulent structures reducing the peak levels of vorticity, 
and 4) increase the levels of numerical dissipation. The dynamic model also shows evidence of providing a 
backscatter of energy into the resolved scales from the sub-grid scales. 

Finally, computing the viscous terms with high-order numerics consistent with the convective terms 
slightly reduces the integrated enstrophy, indicating the numerical dissipation has increased, in comparison 
with using a fourth-order method. 
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